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Abstract 

Background: The molecular circuitry of living organisms performs remarkably robust regulatory tasks, despite the 
often intrinsic variability of its components. A large body of research has in fact highlighted that robustness is 
often a structural property of biological systems. However, there are few systematic methods to mathematically 
model and describe structural robustness. With a few exceptions, numerical studies are often the preferred 
approach to this type of investigation. 

Results: In this paper, we propose a framework to analyze robust stability of equilibria in biological networks. We 
employ Lyapunov and invariant sets theory, focusing on the structure of ordinary differential equation models. 
Without resorting to extensive numerical simulations, often necessary to explore the behavior of a model in its 
parameter space, we provide rigorous proofs of robust stability of known bio-molecular networks. Our results are in 
line with existing literature. 

Conclusions: The impact of our results is twofold: on the one hand, we highlight that classical and simple control 
theory methods are extremely useful to characterize the behavior of biological networks analytically. On the other 
hand, we are able to demonstrate that some biological networks are robust thanks to their structure and some 
qualitative properties of the interactions, regardless of the specific values of their parameters. 



Background 

The complex biochemistry of living organisms very 
often outperforms electrical and mechanical devices in 
terms of adaptability and robustness. Mapping such 
intricate reaction networks to high level design princi- 
ples is the goal of systems biology, and it requires an 
immense collaborative effort among different disciplines, 
such as physics, mathematics and engineering [1]. 

The most classical example of robust molecular circui- 
try is probably given by bacterial chemotaxis [2,3]. The 
action of the flagellar motor of E. coli cells is driven by 
a cascade of signaling proteins, whose active or inactive 
state is determined by the presence of nutrient in the 
environment. Both analysis on a simplified ordinary dif- 
ferential equation (ODE) model [2] and experiments [3] 
showed how the flagellar motion of E. coli presents a 
robustly stable steady state: steps in the nutrient con- 
centration only temporarily alter the motor equilibrium. 
Cells are therefore sensitive to nutrient gradients, but 
always return to their stable motion mode (such 
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property is also referred to as adaptability). Such stable 
steady state is only a function of the concentrations of 
the signaling cascade protein components and a few 
binding rates, and is therefore independent of external 
inputs. Further analysis also demonstrated how integral 
feedback is present in the chemotaxis network, and 
guarantees robustness (perfect adaptation) of the equili- 
brium [4]. 

In this work, we are going to ask a simple question: are 
there biological systems that present structurally stable 
equilibria, and preserve this property robustly with 
respect to their specific parameters? This question has 
been considered before in the literature. For instance, 
through extensive numerical analysis on three-node net- 
works, the authors of [5] have shown that adaptability of 
these systems can be investigated solely based on their 
structure, regardless of the chosen reaction parameters. 
In [6], through numerical exploration of the Jacobian 
eigenvalues for two, three and four node networks, the 
authors isolated a series of interconnections which are 
stable, robustly with respect to the specific parameters. 
Such structures also turned out to be the most frequent 
topologies in existing biological networks databases. 
Numerical simulation has arguably been the most 
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popular tool to investigate robustness of biological net- 
works [7-12]. Analytical approaches to the study of 
robustness have been proposed in specific contexts. A 
series of recent papers [13,14] focused on input/output 
robustness of ODE models for phosphorylation cascades. 
In particular, the theory of chemical reaction networks is 
used in [14] as a powerful tool to demonstrate the prop- 
erty of absolute concentration robustness. Indeed, the so- 
called deficiency theorems [15] are to date some of the 
most general results to establish robust stability of a che- 
mical reaction network. Monotonicity is also a structural 
property that is useful to demonstrate robust dynamic 
behaviors of a class of biological models [16,17]. Robust- 
ness has also been investigated in the context of com- 
partmental models, which are often encountered in 
biology and chemistry [18]. 

In this work, we provide a simple and general theore- 
tical tool kit for the analysis bio-molecular systems. 
Such tools are suitable for the investigation of robust 
stability by means of Lyapunov and set-invariance meth- 
ods. Provided that certain standard properties are veri- 
fied, we demonstrate how a number of well known 
biological networks are asymptotically stable, robustly 
with respect to the model parameters. In some cases, we 
are also able to provide robust bounds on the system 
performance. Our approach does not require numerical 
simulation efforts. The contribution of the paper can be 
summarized as follows. 

♦ The framework we suggest is easy and intuitive for 
biologists to formulate qualitative models without 
the need of exact mathematical expressions and 
parameters. We will propose analytical methods that 
only rely on qualitative interactions between network 
components. 

♦ The properties that can be derived from such 
modeling are, consequently, structurally robust 
because they are not inferred from mathematical for- 
mulas arbitrarily chosen to fit data. 

♦ We suggest techniques based on set-invariance and 
Lyapunov theory, in particular piecewise-linear func- 
tions, to show that such models are amenable for 
robust investigation by engineers and mathemati- 
cians. Such techniques are believed to be quite effec- 
tive and promising in dealing with biological 
robustness [19], [20]. 

♦ We consider several models from the literature, 
reporting the original equations, and rephrasing 
them in our setup as case studies. 

♦ We show how robust certifications can be given to 
important properties (some of which have been 
established based on specific models). 



Methods 

Robustness 

We will consider biological dynamical systems which are 
successfully modeled with ODEs and can be written as: 

x=f(x,u), x(0) = Xq, (1) 

where x is the system state, u models external inputs, 
and both are vectors of appropriate dimensions. Such 
class of models is appropriate for biological systems 
where stochasticity and anisotropy can be neglected. We 
define robustness as follows: 

Definition 1 Let Cbe a class of systems and Vbe a 
property pertaining such a class. Given a family 
T C Cwe say that V is robustly verified by t¥, in short 
robust, if it is satisfied by each element of <¥. 

Countless examples can be brought about families & 
and candidate properties. In this work, we will focus on 
the property of stability, which is an important feature 
for the equilibria of biological networks [1,6,17]. If we 
take into account a linear or linearized dynamical sys- 
tem, we can immediately provide an example that clari- 
fies our definition of robustness [21]. Let C be the class 
of linear differential systems and & the family of second 
order systems described by 
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with positive and constant coefficients a, b, c, d. 
Assume V = asymptotic stability. Then we can say that V 
is robust. The situation is different if we admit that a(t), 
b(t), c(t), d(t) can vary with time, yielding a system 
which is possibly unstable. 

If one is interested in the global system behavior, Lya- 
punov functions are a powerful tool providing sufficient 
conditions for stability. Given an equilibrium point x, 
any convex function V(x — x) > 0 for x =/x and zero at 
the origin is a candidate Lyapunov function. Iff(x, u) is 
continuous, and V (♦) is smooth, then V (♦) is a Lyapu- 
nov function if: 

V(x — x) = VV(x — x)f(x, u) < k(x — x), 

where u is fixed and k (♦) is a negative definite func- 
tion (i.e. k (♦) < 0 on all its domain, except for n (0) = 
0). 

Non-smooth Lyapunov functions 

The concept of Lyapunov derivative can be generalized 
when the function V (♦) is non-smooth. For instance, 
consider the convex function: 

V(x — x) = max — x), i = I, . . . , N, 
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where each V/(-) is smooth and convex, and assume 
that V (•) is positive definite. The set of active functions 
is never empty and is defined as: A = {i : Vi(x — x) = V}. 
If we define the generalized Lyapunov derivative as: 

D + V{x — x) = maxVVi(x — x)f(x), 

ieA 

then the condition for stability becomes: 
D + V(x — x) < k(x — x), /c(-) negative definite. 



Positively invariant sets 

We are interested in cases where the trajectories of sys- 
tem (1) remain trapped in bounded sets at all times, 
therefore behaving consistently with respect to some 
desired criterion. 

We say that a subset <S of the state space is positively 
invariant if x(Q) e S implies that also x(t) e S for all t > 
0. The following theorem (which relies on the concept 
of Lyapunov function) provides a general necessary and 
sufficient condition for a set to be invariant. 

Theorem 1 (Nagumo, 1943) Assume that system (1) 
(for a fixed constant input u) admits a unique solution. 
Consider the set: 

S = [x e R n : s f (x) < cr ir i=l,...,r}, 

where s t are smooth functions, and Ct are given con- 
stants. Assume that Vsj ^ 0, Wx e dS. The set of active 
constraints is l(x) = {i : s;(x) = or f }, and is non-empty 
only on the boundary of S. Then the set Sis positively 
invariant if and only if 

Vsi(x)f{x, u) < 0, Vac £ dS, and i e T{x). 

For instance, if our constraining functions are linear, s T x 
< O", the Nagumo conditions are s T f(x, u) < 0. We refer 
the reader to [22] for further details on positively invariant 
sets; more recent works on this topic are [23] and [24]. 

Structural robustness investigation for biological 
networks 

Let us begin with a simple biological example. Consider 
a protein x lf which represses the production of an RNA 
species x 2 . In turn, x 2 can be the target of another RNA 
species u 2 (and form an inactive complex to be 
degraded) or it can be translated into protein x 3 . A stan- 
dard dynamical model [25] of this process is: 

X\ = U\ — bi\Xi, 

x 2 = d 2 i(xi) - b 22 x 2 - b 2ll2 x 2 u 2 , d 2 i(xi) = 1 , (2) 

1 + x{ 

^3 = #32^2 — ^33^3- 



RNA species x 2 determines the production rate of pro- 
tein x 3 by indexing the corresponding reaction rate as 
a 32 . Following the standard notation in control theory, 
we assume that the production rate of protein Xi is dri- 
ven by some external signal or input u±, and that RNA u 2 
also acts as an external input on RNA x 2 . We assume 
that all the system parameters are positive and bounded 
scalars. Terms a^ are first order production rates: species 
i is produced at a rate which is linear in species y; b ih 
denote in this case first order degradation rates. The 
term ^21(^1) is a well known Hill function term [26]. The 
stability properties of this small network can be immedi- 
ately assessed: x x will converge to its equilibrium 
X\ = U\\b\\. Similarly, x 2 = d 2 i(xi)/(b 22 + b 2u2 u 2 ), 
X3 = ^32^2/^33. Regardless of the specific parameter 
values, and therefore robustly, the system is stable. The 
equilibrium x\ could grow unbounded with u lt however 
x 2 is always bounded. 

We remark that the knowledge of functions a^x, b ih x 
and d{-) is not necessary at all: the previous conclusions 
can be easily derived by the qualitative information that 
d 2 i is strictly decreasing and asymptotically converging 
to 0, while buXi, b 22 x 2l b 2ll2 x 2 u 2 , a 32 x 2 and b 33 x 3 are 
increasing. 

It is appropriate at this point to outline a series of 
general assumptions that will be useful in the following 
analysis. 

We will consider a class of biological network models 
consisting of n first order differential equations 

(3) 

+ ^2c is (x s ) + ^2da(xi), 

seQ leVi 

where x it i - 1,..., n are the dynamic variables. For 
the sake of notation simplicity, we are not denoting 
external inputs with a different symbol. Inputs can be 
easily included as dynamic variables x u = w u {x Ul t) 
which are not affected by other states and have the 
desired dynamics. The sets A\, B\, C\, V\ denote the 
subsets of variables affecting x t . The different terms 
in equation (3) are associated with a specific biologi- 
cal and physical meaning. The terms fji- , ♦) are asso- 
ciated with production rates of reagents, typically, 
these functions are assumed polynomial in their argu- 
ments; similarly, terms g ih (♦ , ♦) model degradation or 
conversion rates and are also likely to be polynomial 
in practical cases. Finally, terms c(-) and d{>) are asso- 
ciated with monotonic nonlinear terms, often given 
by Michaelis-Menten or Hill functions [26]. We 
assume that system (3) satisfies the following 
assumptions: 
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Al (Smoothness) Functions fj (• , •)> g ih (♦ , c is (♦) <2#<i 
d// (♦) ^re unknown nonnegative continuously differenti- 
able functions. 

A2 fj (xu 0) = 0 and g ih (x b 0) = 0, V#. 

A3 Functions fij (x i} xj) and g^ (x i} Xyj), are strictly 
increasing in Xj and Xh respectively. 

A4 (Saturation) Functions c is (x s ) and du(xj) are nonne- 
gative and, respectively, non- deer easing and non-increas- 
ing Moreover c is (°°) > 0 and dn(0) > 0. 

A5 Functions g ih (- , ♦) are null at the lower saturation 
levels : ^(0, x\j) - 0, Va^. 

In view of the nonnegativity assumptions and 
Assumption 5, the general model (4) is a nonlinear posi- 
tive system, according to the next proposition, and its 
investigation will be restricted to the positive orthant. 

Proposition 1 The nonnegative orthant x t > 0 is posi- 
tively invariant for system (4). 

Given the above assumptions, we can write equation (3) 
in an equivalent form. First of all, in view of A1-A3, we 
can write: fj(x if xj) - a(x i} xj)Xp gih(Xp Xp) = b(x i} x^x^, with 

f > fiMuXj) gih{Xi,X h ) 

aij[Xi,Xj) = — — , and Dij(Xi,Xh) = , 

Xj Xfo 

The above expression is always valid due to the 
smoothness assumption Al (see [18], Section 2.1). 

Additionally, assumption A5 requires that b ih (0, Xpj) = 
0, Vxh, for i * h. Once we adopt this notation, we can 
rewrite model (3) as follows: 

Xi{f) = ^ aij(xi,Xj)xj — ^ bih[xi,Xh)xh+ 

+ ^2 c is( x s) + XI faM, i=h2, n. 

seCi leVi 

To simplify the notation, we have considered func- 
tions depending on two variables at most. However, we 
can straightforwardly extend assumptions A1-A5 to 
multivariate functional terms in equation (3). In turn, 
the model structure (4) can be easily generalized to 
include terms as a{x t , Xp xj <} ...), b(x b Xp xp 
x h ..)> d(x if x p x h ...). 

If we restrict our attention to the general class of 
models (4), under assumptions A1-A5, we can proceed 
to successfully analyzing the robust stability properties 
of several biological network examples. 

The structural analysis of system (4) can be greatly 
facilitated whenever it is legitimate to assume that func- 
tions a, b, c d have certain properties. For the reader's 
convenience, a list of possible properties is given below. 
Given a general function^): 



PI f(x) = const > 0 is nonnegative-constant. 

P2 fx) = const > 0 is positive-constant. 

P3 / (x) is sigmoidal: it is non-decreasing, fO) - f '(0) = 
0, if 0 <f(°°) <°° and its derivative has a unique maxi- 
mum point, f(x) < f(x)for some x > 0. 

P4/(#) is complementary sigmoidal: it is non-increas- 
ing, 0 <f(0),f(0) = 0, foo) = 0 and its derivative has a 
unique minimum point. In simple words, f is a CSM 
function ijffO) - fx) is a sigmoidal function. 

P5 f (x) is constant-sigmoidal, the sum of a sigmoid 
and a positive constant. 

P6 / (x) is constant-complementary-sigmoidal, the sum 
of a complementary sigmoid and a constant. 

P7 f(x) is increasing-asymptotically-constant: f(x) > 0, 
0 <f(°°) <°° and its derivative is decreasing. 

P8 / (x) is decreasing-asymptotically-null: f(x) < 0, / 
(oo) = 0 and its derivative is increasing. 

P9 / (x) is decreasing-exactly-null: f(x) < 0, for x <x 
and f (x) - 0 for x > xfor some x > 0. 

P10 f(x) is increasing-asymptotically-unbounded: f(x) 

> 0,/(°o) = + oo. 

As an example, the terms d{>) and c(-) are associated 
with Hill functions, which are sigmoidal and comple- 
mentary sigmoidal functions. A graphical sketch of their 
profile is in Figure 1C and ID. 

Network graphs 

Building a dynamical model for a biological system is 
often a long and challenging process. For instance, to 
reveal dynamic interactions among a pool of genes of 
interest, biologists may need to selectively knockout 
genes, set up micro RNA assays, or integrate fluorescent 
reporters in the genome. The data derived from such 
experiments are often noisy and uncertain, which 
implies that also the estimated model parameters will be 
uncertain. However, qualitative trends can be reliably 
assessed in the dynamic or steady state correlation of 
biological quantities. 

Graphical representations of such qualitative trends 
are often used by biologists, to provide intuition regard- 
ing the network main features. We believe that, indeed, 
such graphs may be useful even to immediately con- 
struct models analogous to (3). We propose a specific 
method to construct such graphs: the biochemical spe- 
cies of the network will be associated to the nodes in 
the graph, the qualitative relationships between the spe- 
cies will be instead associated with different types of 
arcs: in particular, the terms of a, b, c and d will be 
represented as arcs having different end-arrows, as 
shown in Figure 1A. These graphs can be immediately 
constructed, by knowing the correlation trends among 
the species of the network, and aid the construction of a 
dynamical model. For simple networks, this type of 
graph may provide intuition regarding their behavior 
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Figure 1 Graphical representation of biological networks. A. The arcs associated with the functions a, b, c and d. We will use dashed arcs, 
connecting to arcs of the type a and b to highlight that the associated function is nonlinearly dependent on a species of the network: in the 
example above, o 3] = o 3] {x 2 ). B. The graph associated with equations (2); external inputs are represented as orange nodes. C. Examples of 
sigmoidal functions. D. Examples of complementary sigmoidal functions. In our general model (4), functions d{) and c(-) are naturally associated 
with Hill function terms. 

V ) 



and may facilitate their structural robustness analysis. 
For instance, the graph associated to equations (2) is 
shown in Figure IB. Throughout the paper, we will con- 
sider similar case studies and use their graph representa- 
tion as a visual support for the analysis. 

Remark 1 In this work, properties such as positivity, 
monotonicity, boundedness and other functional charac- 
teristics are labeled as "qualitative and structural prop- 
erties"^]. Through such properties, we can draw 
conclusions on the dynamic behaviors of the considered 
systems without requiring specific knowledge of para- 
meters and without numerical simulations. However, it 
is clear that our approach requires more information 
than other methods, such as boolean networks and other 
graph-based frameworks. 

Investigation method 

The main objective of this work is to show that, at least 
for reasonably simple networks, structural robust stabi- 
lity can be investigated with simple analytical methods, 
without the need for extensive numerical analysis. We 



will suggest a two stage approach: 

♦ Preliminary screening: establish essential informa- 
tion on the network structure, recognizing which 
properties (such as P1-P10) pertain to each link. 

♦ Analytical investigation: infer robustness properties 
based on dynamical systems tools such as Lyapunov 
theory, set invariance and linearization. 

Results and Discussion 

In this section we will analyze five biological networks 
as case studies. Three of such examples, the L-arabinose, 
the sRNA and the Lac Operon networks, model the 
interaction and control of expression of a set of genes. 
The cAMP and the MAPK pathways are instead signal- 
ing networks, namely they represent sets of chemical 
species interacting for transmission and processing of 
upstream input signals. These networks are all well- 
known in the literature, and have been characterized 
mainly through experimental and numerical methods, 
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although the MAPK pathway, for instance, has been 
thoroughly analyzed using the theory of monotone sys- 
tems [17]. 

We will provide rigorous proofs that these networks 
are either mono or multi-stable in a robust manner. 
Such demonstrations rely on Lyapunov functions and 
invariant sets theory, according to our proposed metho- 
dology. In some cases, we are also able to provide 
bounds on their speed of convergence. 

The L-arabinose network 

The arabinose network is analyzed in [28] as an example 
of feedforward loop. Two genes araBAD and araFGH 
are regulated by two transcription factors, AraC and 
CRP. AraC is a repressor, but turns into an activator 
when bound to the sugar L-arabinose. CRP is an activa- 
tor when bound to the inducer cyclic AMP (cAMP), 
which is produced when cells are starving upon glucose 
(not produced during growth on glucose). CRP also 
binds to the araC promoter and enhances transcription 
of AraC, which has a significant basal rate of expression 
(i.e. it is produced by the cell also in absence of inducer 
CRP). A very simple model for this network can be 
derived by defining the state variables X\ and x 2 , respec- 
tively the concentrations of the transcription factor 
AraC and of the output protein araBAD. The concentra- 
tion of the transcription factor CRP is considered an 
external input u: 

xi=pi+ Pif{u,K UXl ) - aiXi, 

x 2 = Pif{u,K UX2 ) -f(xi,K XlX2 ) - a 2 x 2 , 

where a x , a 2 are the degradation and dilution rates of 
Xi, x 2 respectively. The basal production rate of X\ 
(AraC) is pi. The activation pathways are modeled by 
Hill functions f (u, K) = u H l(K H + u H ), where H is the 
Hill coefficient and Kg are the activation thresholds. The 
model can be recast into the general structure (4), 
which includes model (5) as special case:: 

xi = a +ciu(u) - b n x lf 

( \ u (6) 
x 2 = C 2u i[U, X\) — U 22 X 2 , 

where u is nonnegative- constant, c lf b X \ and b 22 are 
positive-constant, while Ci u (u) and c 2u i(u) are sigmoidal 
with respect to u, the latter increasing with respect to x ± . 
The graph representation of this network is in Figure 3A. 

For this elementary network the analysis is straightfor- 
ward. Variable X\ is not affected by x 2 . Since Ci u (u) is 
bounded, X\ is also bounded and converges to an equili- 
brium point X\ (u) which is monotonically increasing in 
u. In turn, x 2 is also positive and bounded for any value 
of u and stably converges to a unique equilibrium point 
x 2 , which is a monotonically increasing function of u 
(partially activated by X\(u)). The positive term c x 



prevents Xi(t) and x 2 (t) from staying at zero. It is worth 
remarking that the hierarchical structure of this network 
greatly facilitates the analysis; equilibria can in fact be 
iteratively found and their stability properties 
characterized. 

The sRNA pathway 

Small regulatory RNAs (sRNA) play a fundamental role 
in the stress response of many bacteria and eukaryotes. 
In short, when the organism is subject to a stimulus 
that threatens the cell survival, certain sRNA species are 
transcribed and can down-regulate the expression of 
several other genes. For example, when E. coli cells are 
lacking a source of iron, the sRNA RyhB (normally 
repressed by the ferric uptake regulator Fur) is 
expressed and rapidly induces the degradation of at least 
other 18 RNA species encoding for non-essential pro- 
teins which use up Fe molecules. This allows essential 
iron-dependent pathways to use the low amount of 
available iron. Quantitative studies of the sRNA path- 
ways have been carried out in [29-31]. Let us define X\ 
as the RNA concentration of the species which is tar- 
geted by the sRNA and x 2 as the concentration of 
sRNA. The model often proposed in the literature is: 

k\ = a\ — fiiXi — kx\X 2 , 
X 2 = ot 2 — filX 2 — kX\X 2 , 

where a lt a 2 are the transcription rates of x x and x 2 
respectively, p lf p 2 are their degradation rates (turnover), 
and h is the binding rate of the species X\ and x 2 . The 
formation of the inactive complex X\ ' x 2 corresponds to 
a depletion of both free molecules of x x and x 2 . If a x 
<a 2 the pathway successfully suppresses the expression 
of the non-essential gene encoded by X\. This model 
can be embedded in the general family: 

X\ — C\ — ^n^i — b\ 2 {X\,X 2 }X 2 , 

x 2 = c 2 — b 22 x 2 — b 2 i(xi,x 2 )xi, 

by setting b X2 - kxi and b 2X - kx 2 (note that ^12(0) = 
&2i(0) as required). From our list of properties: c lf c 2) 
bn and b 22 are positive- constant; bi 2 {xi, x 2 ) and b 2 i (xi, 
x 2 ) are increasing-asymptotically-unbounded in both 
variables; and ^12(^1, x 2 )x 2 = b 2 \{x 1) x 2 )xi at all times. 
This network can be represented with the graph in 
Figure 2A. 

The sRNA system is positive, because the nonnegativ- 
ity Assumptions 1 and 4 are satisfied. The preliminary 
screening of this network tells us that each variable pro- 
duce an inhibition control on the other, which increases 
with the variable itself. In other words %\ is "less toler- 
ant" to an increase of x 2 if the latter is present in a large 
amount. This means that the sum x x + x 2 is strongly 
kept under control. Also the mismatch between the two 
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Figure 2 The sRNA network. A. The graph associated with the sRNA network B. Sectors, Lyapunov function level curves (orange) and 
qualitative behavior of the trajectories (green) for the sRNA system 



variables is controlled. 1 To prove stability of the 
(unique) equilibrium x, we will use the 1-norm as Lya- 
punov function V(x — x) = \\x — x\\i (see Figure 2B). This 
choice has a remarkable interpretation: denoting by 
E = (x\ — X\ ) + (x 2 — x 2 ) and A = (x\ — X\ ) — (x 2 — x 2 ) 
the sum and the difference of the two variables (referred 
to the equilibrium) we have 



V(x — x) 



■ x\ 



\X\ — X\ I + \x 2 
max{A, £}, 



■x 2 \ 



thus the function represents the worst case between 
the sum and the mismatch. 




Figure 3 Graphs associated with case studies. A. The graph associated with the L-arabinose network, external inputs are represented as 
orange nodes. B. The graph associated with the cAMP pathway. C. The graph associated with the lac Operon network. D. The graph associated 
with the MAPK signaling pathway. 
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The following proposition shows that the sRNA 
pathway is a typical system in which robustness is 
structurally assured. We report the full demonstration 
of this proposition, because its steps and the techni- 
ques used are a model for the subsequent proofs in 
this paper. 

Proposition 2 The variables of system (8) are 
bounded for any initial condition #i(0), #2(0) > 0. The 
system admits a unique asymptotically stable equili- 
brium point x = (xi,x 2 ) T and the convergence is expo- 
nential: 

\\x{t) - x\\ l < e'P* || jc(0) -x\\ v (9) 

for some /3 > 0 and any Xi(0) > 0, x 2 (0) > 0. Moreover, 
no oscillations are possible around the equilibrium, in 
the sense that the condition X\(f) = X\or x 2 {t) = x 2 occurs 
at most once. 

Proof To prove boundedness of the variables we need 
to show the existence of an invariant set 

S = [x\ > 0, x 2 > 0 : X\ + x 2 < k}. 

Proposition 1 guarantees that the positivity constraints 
are respected. Then we just need to show that the con- 
straint Xi + x 2 ^ k cannot be violated for sufficiently 
large k > 0. The derivative of function s (x lf x 2 ) = x x + 
x 2 is 

s(Xi, X2) = X\ + X2 

= C\ — bi\Xi — fei2(Xi,X2)^2 + C2 — ^22x2 ~ ^21 (^1/^2^1 
< C\ — b\\X\ + C2 — \)22^2 

<C\+C2 — min{bn, ^22}(^i + ^2) 
= c\ + C2 — m'm{bii,b22}s{xi,X2). 

Define k, = (c ± + c 2 )/min {bn, b 22 } then for s(x lf x 2 ) >k 
the derivative becomes negative so s(x lf x 2 ) cannot 
exceed k (See Theorem 1). 

Boundedness of the solution inside a compact set 
assures the existence of an equilibrium point. Let 
(xi,x 2 ) be any point in which the following equilibrium 
conditions holds: 

c\ — b\\X\ — ^12(^1/^2)^2 = 0, 
c 2 — b 22 x 2 — ^21(^1/^2)^1 = 0. 

The behavior of the candidate Lyapunov function: 

V(xi,x 2 ) = \x\ — X\ \ + \x 2 — x 2 \ 

= max{=b(xi — Xi) ± {x 2 — x 2 )}, 

will be examined in the different sectors represented 
in Figure 2B. Let us start by considering the sector 
defined by X\ > X\ and x 2 > x 2 (APB in Figure 2B) for 
which V(xi,x 2 ) = [x] — X]) + (x 2 — x 2 \ In such a sector 
the Lyapunov derivative is: 



D + V(xi,x 2 ) = [1 1] [*] 

= C\ + C 2 — ^11^1 — b 22 X 2 — 

— ^12(^1/^2)^1 — ^21 (^-1 / ^2)^2 

= -b n (xi - xi) - b 22 [x 2 - x 2 )- 

— [^12(^1/^2)^1 — ^12(^1/^2)^1] 

— [^21 C^-i/ ^2)^2 — ^21(^1/^2)^2] 

< — ^11(^1 — ^1) — b 22 {x 2 — ^2)/ 

where we have subtracted the null terms (10) 
and where we have exploited the fact that b 12 {xi, x^x^ 
= b 2 i{xi, x 2 )x 2 is increasing in both variables. 
The inequality (CPD in Figure 2B) 
D + V(xi,x 2 ) < bn(xi —Xi) + b 22 [x 2 — x 2 ) < 0 can be 
similarly proved to hold in the sector X\ < X\ and 
x 2 < x 2 . 

Consider the sector defined by X\ > X\ and x 2 < x 2 
(DPA in Figure 2B) for which V(xi,x 2 ) = X\ — X\ in such 
a sector the Lyapunov derivative is: 

D + V(xi,x 2 ) = [1 -i][y 

= C\ — c 2 — b\\X\ + b 22 x 2 — 
— ^12(^1/^2)^1 + b 2 \{x 2l x 2 ^x 2 

= 0 by assumption 

= —bn (xi — X\) + b 22 (x 2 — x 2 ) < 0. 

Note that in the last step we have added and sub- 
tracted the null terms (10). In the opposite sector (BPC 
in Figure 2B) X\ < X\ and x 2 > x 2i we can prove that 
D + V(xi,x 2 ) = +bn(xi — Xi) — b 22 {x 2 — x 2 ) < 0. We just 
proved that 

D + V(xi,x 2 ) < —[bn \x\ — X\ \ + b 22 \x 2 —x 2 \] 

< -fi[\xi -xi| + \x 2 -x 2 \] 

< -pv(xi,x 2 ), 

with /3 - min{Z?n, ^22}- This implies (9) and the 
uniqueness of the equilibrium point. 

We finally need to show that there are no oscillations. 
To this aim we notice that the sectors DPA, X\ > X\ and 
x 2 < x 2 , and its opposite 

CPB, Xi < X\ and x 2 > x 2 , are both positive invariant 
sets. 

We can apply Nagumo's theorem: consider the half- 
line PA originating in P, where x 2 = x 2 and Xi > X\. 
Therefore we have that (again by adding the null term 
in (10)): 

x 2 — c 2 — b 22 x 2 — b 2 \ {x\, x^Xx — 
— c 2 — b 22 x 2 — ^21(^1/^2)^1 
= —[^21(^1/^2)^1 — ^21(^1/^2)^1] ^ 0. 
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Similarly, on half-line PD where X\ = X\ and x 2 < x 2 , 
by considering (10) we derive 

X\ = — b\ 2 {x\,x 2 }x 2 + b\ 2 {x\,x 2 }x 2 ^ 0 

hence the claimed invariance of sector DPA. The 
proof of the invariance of sector CPB is identical 

Remark 2 Note that the constructed Lyapunov function 
\\x — x\\\does not depend on the system parameters. This 
fact can be used to prove that if the transcription rates c x 
(t) and c 2 (t) are time-varying, but bounded, we have con- 
vergence to a neighborhood whose amplitude, obviously, 
depends on the bounds of c x (t) and c 2 (t). It is realistic to 
assume that the transcription rates vary over time: for 
instance, if the environmental conditions change, the cell 
may need to down or up-regulate entire groups of tran- 
scripts and therefore increase or decrease c 2 (t). 

The following corollary evidences the positive influence 
of c 2 , which is positive over x 2 and negative over %\. 

Corollary 1 Assume that #i(0), x 2 (0) is at the steady 
state corresponding to i\and c 2 . Consider the new input 
d = c\(keeping c\ = c-J. Then the system converges to a 
new equilibrium with X\ < X\{0)and x 2 > x 2 (0). There is 
no undershoot, respectively, overshoot. 

Proof The steady state values X\ and x 2 are respec- 
tively monotonically decreasing and increasing functions 
of c 2 . Indeed, consider the steady-state condition 

C\ = b\X\ + bi 2 {xi,x 2 )x 2l 
c 2 = b 2 x 2 + b 2 \ {x,\, x^xi, 

and regard it as a differentiable map (x\, x 2 ) — > {c\ f c 2 ). 
By the uniqueness proved in Proposition 2 the map is 
invertible. The Jacobian of the inverse map is the inverse 
of the Jacobian 



9(^12^2) 3(^12^2) 



dX\ 
d(b2i*i) 
dX\ 



D 2 + 



dx 2 

3 (621*1 ) 



det(J) 



b 2 + 



dx 2 

d(b 2 ixi) d{b i2 x 2 ) 



dx 2 
d{b 2 iXi) 

dX\ 



dx 2 

d(bi 2 x 2 ) 



dx 



1 - 1 



where 



a +m u u u a (^i2^2) u d{b 2 iXi) 
det(J) = b\b 2 + b 2 — + b\ — > 0 



dX\ dx 2 

(keep in mind that b 2 i(xi, x 2 )xi = bi 2 (xi, ^2)^2)- The 
sign of the entries in the second column are negative 
and positive respectively, therefore, the steady-state 
values x\ and x 2 are decreasing and increasing functions 
of c 2 . 

The absence of overshoot and undershoot is an 
immediate consequence of the invariance of the sector 
X\ > X\ and x 2 < x 2 previously proved. 



Obviously, decreasing c 2 increases x x and decreases x 2 
and the same holds if we commute 1 and 2. It is worth 
noting that the same conclusions regarding the lack of 
multistability and oscillations for the sRNA pathway 
may be reached by qualitative analysis of the system's 
nullclines. 

The cAMP dependent pathway 

The cyclic adenosine monophosphate (cAMP) pathway 
can activate enzymes and regulate gene expression 
based on sensing of extracellular signals. Such signals 
are sensed by the G protein-coupled receptors on the 
cell membrane. When a receptor is activated by its 
extracellular ligand, a series of conformational changes 
are induced in the receptor and in the attached intracel- 
lular G protein complex; the latter activates adenylyl 
cyclase, which catalyzes the conversion of ATP in 
cAMP. In yeast, cAMP causes the activation of the pro- 
tein kinase A (PKA), which in turn regulates the cell 
growth, metabolism and stress response. Stochastic 
models are usually proposed for numerical analysis of 
the cAMP pathway. However, the cAMP pathway com- 
ponents in yeast are present in high numbers and a 
deterministic modeling approach is adopted in [31]. In 
such paper, the pathway is decomposed in several mod- 
ules, here we consider the simplified cAMP Model A, 
which focuses on the parts of the pathway best charac- 
terized in the literature: 



X\ = kp^f — X\)U — UrXi, 



X 2 '■ 



kp^x^ — ^2)^3 — kj^X 2 , 
Of 3 + k A Xi V maxPl x 2 x 3 



(ii) 



1 + kjX 2 K]\4p x + -KmP 2 + x 3 



where X\ is the concentration of active G protein, x 2 is 
the concentration of active PKA protein, x 3 is the con- 
centration of cAMP and u is the concentration of glu- 
cose input to the network. The parameters V max p 1 and 
VmaxP 2 model the "feedback" effect introduced by two 
phosphodiesterases (Pdelp and Pde2p) on the cAMP 
concentration. The numerical analysis in [31] typically 
shows that the cAMP concentration (x 3 ) responds with 
a large overshoot to steps in the glucose (u, input) con- 
centration. We will analytically explore the dynamic 
behavior of x 3 , showing that under certain assumptions, 
a bounded overshoot is a robust characteristic in the 
system. The parameters h F and h R are forward and 
reverse reaction rates for the formation of active X\ and 
x 2 . Mass conservation allows to express the active pro- 
tein amounts as a function of the total number of mole- 
cules, xf l = x™ actwe + x^ The nonlinear expressions in 
equation x 3 are derived by Michaelis-Menten enzymatic 
steps. We can re-write the above equations according to 
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the general model (4): 
k\ = ai u (xi)u — b\\X\, 

%2 = ^23(^2)^3 — ^22^2/ (12) 
^3 = d 32 (x 2 ) + ^3l(^2)^l — ^32fe)^2 — b 33 (x 3 )x 3l 

where u is the external signal and where b 23 = 0 for x 2 
= 0 and b 32 = 0 for x 3 = 0. A qualitative graphical repre- 
sentation of this network is in Figure 3B. 

Our preliminary analysis allows us to assume: a Xw a 23 : 
deer easing- exactly -null with threshold values x l ° l and 
#2*' ^32> ^3i : deer easing-asymptotically -null, b 32 and g 33 
- ^33fe)^3 : increasing-asymptotically-constant; b 11} b 22 
are positive-constant. 

It is immediate to notice that for constant u, X\ 
robustly converges asymptotically to its equilibrium 
value such that 



The solution X\ = %{u) of the previous equation is 
uniquely defined for each u since the function ^(xi) on 
the right is strictly increasing and grows to infinity, pre- 
cisely lim Xl ^ x tot£ -1 (xi) = +00. Biologically, this means 
that external glucose signals are mapped to internal 
active G-protein concentration with a bijection, before 
saturating. 

Also we have to note that the model is consistent with 
mass conservation: since ai u (xi) and a 23 (x 2 ) are zero 
above the thresholds and x\° l , we have that x\ < 0 
and x 2 < 0 for x\ > and x 2 > x l 2 , respectively; there- 
fore we assume X\(f) < x l [ l , x 2 {t) < x l 2 , for all t > 0. 

Proposition 3 There exists an equilibrium for system 
(12) if and only if 

d 31 {xf) + a 3 i{xf)x x < lim [b 32 (x 3 )xf + ^33(^3)^3]/ (13) 

where X\ = £ (u)as previously defined . All the equili- 
brium values X\ = § {u\ x 2 and x 3 are increasing functions 
of u. If condition (13) is satisfied , the equilibrium is 
unique and locally stable. 

The previous proposition assures only local stability, 
but this result can be extended to global stability. To 
this aim, we will assume that X\ is at its equilibrium 
value X\. Furthermore, under a suitable condition a per- 
formance bound on the transient values of x 3 (f) can be 
given. 

Proposition 4 Assume that x x has reached its steady 
state x\. Then, the unique equilibrium point is globally 
attractive for any initial condition x 2 (0), x 3 {0) > 0. More- 
over, assume that 

l 3 = lim b 3 3(x 3 )x 3 > d 32 (0) + a 3 i(0)£, (14) 



then we can give the following bound for the transient 
ofx 3 {t) 

x 3 (t) < max{x 3 (0),d 3 2(0) + a 3 i(0)^}. (15) 

The proof can be found in Section S{II of the Addi- 
tional File 1. 

Remark 3 The condition (14) has the following inter- 
pretation. It basically states that the inhibiting term b 33 
{x 3 )x 3 at "full force" (x 3 suitably large) dominates the 
activating term d 32 (x 2 ) + a 3X [x 2 )^ when x 2 is small. Note 
that, indeed, the feedback terms modulated by the two 
phosphodiesterases act in a complementary manner, in 
order to maintain a bounded concentration of cAMP in 
the cell. 

Remark 4 The system, even if initialized with small 
values x 2 (0) and x 3 (0), may exhibit a spike of cAMP x 3 
which is bounded by (15), if condition (14) is satisfied. If 
x 3 (0) is small, then the bound is d 32 (0) + <2 31 (0)<f (u): the 
amplitude of the spike is, in general, an increasing func- 
tion of the glucose concentration u. If condition (14) fails, 
then (see Figure S2 in the Additional File) the spike of x 3 
(t) can be arbitrarily large; thus condition (14) can be 
seen as a threshold. 

The Lac operon 

This genetic network was originally studied by Monod 
and Jacob [33] . The natural nutrient for E. coli bacterial 
cells is glucose, which is metabolized by enzymes nor- 
mally produced by the bacteria. When glucose is absent, 
but the allolactose inducer is present in their environ- 
ment, E. coli activates a set of genes that will regulate 
the lactose intake and breakdown. In particular, the cells 
start producing a permease protein, which binds to the 
cell membrane and increases the inflow of lactose; and 
cells also start producing the /3-galactosidase protein, 
which converts lactose in allolactose. 

In this paper we will consider the deterministic model 
proposed in [34]. This simple model does not capture 
the stochasticity of this genetic circuit, but it does 
explain the bimodal behavior of the system. Such beha- 
vior is observable experimentally: within the same popu- 
lation, the operon can be either induced or uninduced. 
Our analysis shows that for low or high intracellular 
inducer concentrations, the system is monostable and 
respectively reaches an uninduced or induced equili- 
brium; however, at intermediate inducer concentrations 
the system becomes multi-stable. 

The state variables of the ODE model we will study 
are the concentration of nonfunctional permease protein 
Xi, the concentration of functional permease protein x 2 ; 
the concentration of inducer (allolactose) inside the cell 
x 3 , and the concentration of /3-galactosidase # 4 , a quan- 
tity that can be experimentally measured. The 
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concentration of inducer external to the cell is here 
denoted as an input function u. 

X\ = /lfe) — 8\Xi, 

x 2 = B\X\ — 8 2 x 2 , 

H (16) 

X3 = \fl{u) — h{x^) ]X 2 + fi 2 U — 83X3, 

X4 = yf i(x3) — 84X4., 

where /5 lf fi 2 , S lf 5 2 , S 3 and / are constants and^ are func- 
tions that are experimentally measurable. In particular, at 
low inducer concentrations, j x ~ k\ + k 2 x 3 + ^3X3, where 
/c/s are constant; at high x 3 concentrations f saturates. The 
functions^ and ^3 are assumed to depend hyperbolically on 
their arguments. According to the proposed setup, the pre- 
vious equations can rewritten as follows: 

X\ — ^13(^3) — h\\X\, 
x 2 = d 2 \X\ — b 22 x 2 , 

(17) 

X3 = Cl3 2 (u)x 2 — b3 2 (X3 )X 2 + C3 U U — ^33X3, 
X4 = ^43(^X3) — ^44X4, 

where c 13 {x 3 ) = fi(x 3 ), b ±1 = S 1} a 2l = p lf b 22 = S 2 , a 32 
(u) =f 2 {u) =, b 32 {x 3 ) =f 3 {x 3 ), c 3u = fi 2 , b 33 = 5 3 , c 43 {x 3 ) = 
7fi(x 3 ) and & 44 = S 4 . This corresponds to the network in 
Figure 3C. 

From our preliminary analysis step: c 13 is constant-sig- 
moidal, a 32 {u) and b 32 {x 3 ) are increasing-asymptotically- 
constant, and the remaining functions a 21 , b llf b 22 and 
b 33 are positive-constant. 

We can start to study this network without any speci- 
fic knowledge of the parameters in equations (17). First 
of all, as evident in Figure 3C, note that the /3-galactosi- 
dase concentration x 4 does not affect any other chemical 
species: therefore, the fourth equation can be considered 
separately. As long as the inducer concentration of x 3 
within the cell reaches an equilibrium X3, X4 converges 
to x 4 = C43(x 3 )/?744. Therefore, we can restrict our atten- 
tion to the first three equations; this is consistent with 
the model proposed in [35,36]. From now on we will 
consider this reduced model (see Section S-III of the 
Additional File), neglecting the linear term c 3u u as in 
[35,36]. We will not introduce delays in our model, as 
done in [37]. Our preliminary screening also shows that 
the evolution of this system is necessarily bounded. 
Indeed Xi receives a bounded signal from x 3 and the 
degradation term -bx\X\ keeps X\ bounded. In turn, x 2 
remains bounded. The inducer concentration x 3 receives 
a bounded signal form u and x 2 ; therefore x 3 stays 
bounded as well, being both a 32 (u) and ^32(^3) bounded. 

The following proposition evidences that fundamental 
results can be established starting from our general fra- 
mework. These results are consistent with the findings 
in [36], whose analysis relies on assuming Hill-type 
functions in the model. 



Proposition 5 For any functional terms in Equations 
17, satisfying the general assumptions formulated above, 
the system admits a unique equilibrium for large u > 0 
or small u > 0. 

For some chioces of such functional terms, the system 
may have multiple positive equilibria x^, x B , x c e IR 3 
(typically three) for intermediate values of u. If multiple 
equilibria exist, then they are ordered in the sense that 
< x B < x c ... where the inequality has to be considered 
componentwise. If the equilibria are all distinct, then 
they are alternatively stable and unstable. In the case of 
three equilibria, x*, x B , x c they are stable, unstable and 
stable, respectively. Finally, given any equilibrium point, 
the positive and negative cones x < x* and x > x* are 
positively invariant. 

The proof is given in Section S-III of the Additional 
File. The cone-invariance property implies that the state 
variables cannot exhibit oscillations around their equili- 
bria. For instance, if x A is the first (hence stable) equili- 
brium, given any initial condition upper bounded by x A 
(x(0) x A ) in the domain of attraction, the convergence to 
x A has no overshoot (and if #(0) > x A there is no 
undershoot). 

Remark 5 It is interesting to notice that, due to the 
competition between terms a 32 and b 32 , the considered 
Lac Operon model is not a monotone system according 
to the definition in [16], where a different model was 
considered. 

MAPK signaling pathway 

Mitogen-activated protein (MAP) kinases are proteins 
that respond to the binding of growth factors to cell 
surface receptors. The pathway consists of three 
enzymes, MAP kinase, MAP kinase kinase (MAP2K) 
and MAP kinase kinase kinase (MAP3K) that are acti- 
vated in series. By activation or phosphorylation, we 
mean the addition of a phosphate group to the target 
protein. Extracellular signals can activate MAP3K, which 
in turn phosphorylates MAP2K at two different sites; in 
the last round, MAP2K phosphorylates MAPK at two 
different sites. The MAP kinase signaling cascade can 
transduce a variety of growth factor signals, and has 
been evolutionary conserved from yeast to mammals. 

Several experimental studies have highlighted the pre- 
sence of feedback loops in this pathway, which result in 
different dynamic properties. This work will focus on a 
specific positive-feedback topology, where doubly-phos- 
phorylated MAPK has an activation effect on MAP3K. 
Such positive feedback has been extensively studied in 
the literature, since the biochemical analysis of Huang 
and Ferrell [37,38] on the MAPK cascade found in 
Xenopus oocytes. In this type of cells, Mos (MAP3K) 
can activate MEK (MAP2K) through phosphorylation of 
two residues (converting unphosphorylated MEK to 
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monophosphorylated MEK-P and then bisphosphory- 
lated MEK-PP). Active MEK then phosphorylates p42 
(MAPK) at two residues. Active p42 can then promote 
Mos synthesis, completing the closed positive-feedback 
loop. 

The presence of such positive feedback in the MAPK 
cascade has been linked to a bistable behavior: the 
switch between two stable equilibria in Xenopus oocytes 
denotes the transition between the immature and 
mature state. A standard ODE model for the cascade is 
proposed in [17], where the authors demonstrate bi-sta- 
bility of the system by applying the general theory of 
monotone systems. We adopt such model, which is 
reported below: 
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(18) 



K 8 + z 2 K w + z 2 ' 



where x is concentration of Mos (MAP3K), y x is the 
concentration of unphosphorylated MEK (MAP2K), y 2 is 
the concentration of phosphorilated MEK-P, y 3 is the 
concentration of MEK-PP, z lf z 2 and z 3 are respectively 
the concentrations of unphosphorylated, phosphorylated 
and doubly-phosphorylated p42 (MAPK). Finally, u is 
the input to the system. 

While bi-stability may occur due to other phenomena, 
such as multisite phosphorylation [39], rather than due 
to feedback loops, a large body of literature focuses on 
bi-stability induced by the positive-feedback in the 
Huang-Ferrel model in Xenopus [40,41] reported above. 
In [37] the feedback/ (u) was characterized, through in 
vitro studies, as an activating Hill-function with high 
cooperativity. In [17] instead, f (u) was assumed to be a 
first order linear term in the concentration of active 
MAP3K, x 7 . In Proposition 6, we will explore the effects 
of different qualitative functional assumptions on the 
feedback loop dynamics f(u). We will highlight that the 
system loses its well-known bi-stability not only in the 
absence of feedback, but also when the feedback 
becomes unbounded. An unbounded positive feedback 
would be caused, for instance, by an autocatalytic 



process of MAP3K activation, mediated by active 
MAPK. We choose to rewrite the above model as fol- 
lows: 

%i = a\ 7 [x\)iix 7 + C\o — bn(xi)xi, 
%2 = c 23 (x 3 ) — b 2 i(x 2 )xi, 

%3 = ^3l(^2)^l + £34(^4) — ^31(^3)^1 — ^33(^3)^3/ 

x 4 = a 4 i{x 3 )x\ — b 44 {x 4 )x 4l (19) 
x 5 = c 56 (x 6 ) - b 54 {x 5 )x 4 , 

Xfr = #64(^5)^4 + £67 0*7) — ^64(^6)^4 — ^66(^6)^6/ 

x 7 = a 74 (x 6 )x 4 — b 77 (x 7 )x 7 . 

The term ftx 7 introduces the positive feedback loop 
and represents a key parameter for the analysis to fol- 
low. A preliminary screening of the system immediately 
highlights the following properties. Function bn(xi)xi, 
functions c 23 (x 3 ), b 21 (x 2 ), a 4X (x 3 ) and b 44 (x 4 )x 4 , functions 
c 56(^6)> b 54 {x 5 ), ^74(^6) an d ^77(^7)^7 are increasing- 
asymptotically -constant. Moreover, a 3 i (x 2 ) = b 2 i(x 2 ), c 34 

(x 4 ) = b 44 (x 4 )x 4 , b 3 i(x 3 ) = ^4lfe)> ^33fe)^3 = c 23fe) an< i 

^64(^5) ~ b 54 (x 5 ), c 67 (x 7 ) = b 77 (x 7 )x 7 , b 64 (x 6 ) = a 74 (x 6 ), 
b^ x 6) x 6 - c 56(^6)« We assume c 10 to be a positive- 
constant. 

The graph in Figure 3D can be partitioned considering 
three aggregates of variables, precisely {xi\, £234 = {x 2) 
x 3) x 4 ) and S 567 = {x 5 , x 6 , x 7 }. Signal X\ is the only input 
for S234, signal x 4 is the only input for E 567 . Then x 7 is 
fed back to the first subsystems by arc a 17 . Without the 
positive feedback loop, we will demonstrate that the sys- 
tem is a pure stable cascade. Note also that £234 and 
£ 567 can be reduced since x 2 + x 3 + x 4 = 0 and 
Xs + x 6 + x 7 = 0 and therefore the following sums are 
constant 



x 2 (t) + x 3 (t) + x 4 {t) = k, 
X5(t) + x 6 (t) + x 7 (t) = h, 



(20) 



with h - x 2 (0) + #3(0) + #4(0) and h - x 5 (0) + x 6 (0) + 
x 7 (0). Since x t > 0, all the variables but X\ are bounded. 
The system can be studied by removing variables x 3 = h 
- x 2 - x 4 and x 6 = h - x 5 - x 7 . We must assume that 
Ciolim^^oo^n^ijxi otherwise no equilibrium is possi- 
ble. The following result is proved in Section S-IV of 
the Additional File. 

Proposition 6 For ft = 0 the system admits a unique 
globally asymptotically stable equilibrium. 

For ft > 0, the system may have multiple equilibria, for 
specific choices of the involved functions a, b, c. 

For ft > 0 suitably large and a 17 (xi) lower bounded by 
a positive number, then the system has no equilibria. 

For ft > 0 suitably bounded and ai 7 (xi) increasing, or 
non-decreasing, and bounded, if multiple simple 2 equili- 
bria exist, then such equilibria are alternatively stable 
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and unstable. In the special case of three equilibria, then 
the system is bistable. 

For ft > 0 suitably bounded and a 17 (x\) increasing 
asymptotically unbounded, then the number of equilibria 
is necessarily even (typically 0 or 2). Moreover, if we 
assume that there exists ft* > 0 such that the system 
admits two distinct equilibria for any 0 <ft < ft*, then 
one is stable, while the other is unstable. 

The proof of this last proposition also shows that mul- 
tiple equilibria x A , x B have a partial order: 
x^ <x^ <x^ . . . f x^<x^<x^...x^<x^<Xy... while 
x 2 and X5 have the reverse order x^ >x^ >x 2 . . . and 

x^ ^5 • • • 

Remark 6 The simplest case of constant a 17 has been 
fully developed in [17] 3 and [16], and it turns out that 
the system may exhibit bi-stability for suitable values of 
the feedback strength ft. Here we show that, for constant 
a 17 , bi-stability is actually a robust property. Our results 
are consistent with the fact that the MAPK cascade is a 
monotone system and some of them could be demon- 
strated with the same tools used in [16,17]. With respect 
to such literature, our contribution is that of inferring 
properties such as number of equilibria and mono or bi- 
stability starting from qualitative assumptions on the 
dynamics of the model, without invoking monotonicity. 

Remark 7 Finally, it is necessary to remark that our 
results on the MAPK pathway robust behaviors hold true 
given the model (19) and its structure. Other work in the 
literature shows that feedback loops are not required to 
achieve a bistable behavior in the MAPK cascade [38], 
when the dual phosphorylation and de-phosphorylation 
cycles are non-processive (i.e. sites can be phosphory- 
lated/ de-phosphorylation independently) and distributed 
(i.e. the enzyme responsible for phosphorylation/ de-phos- 
phorylation is competitively used in the two steps). 

Conclusions 

A property is structurally robust if it is satisfied by a 
class of systems of a given structure, regardless the 
choice of specific expressions adopted and of the para- 
meter values in the model. We have considered five 
relevant biological examples and proposed to capture 
their dynamics with parameter-free, qualitative models. 
We have shown that specific robust properties of such 
models can be assessed by means of solid theoretical 
tools based on Lyapunov methods, set-invariance the- 
ory and matrix theory. Robustness is often tested 
through simulations, at the price of exhaustive cam- 
paigns of numerical trials and, more importantly, with 
no theoretical guarantee of robustness. We are far 
from claiming that numerical simulation is useless. It 
it important, for instance, to falsify "robustness conjec- 
tures" by finding suitable numerical counterexamples. 



Furthermore, for very complex systems in which analy- 
tic tools can fail, simulation appears be the last resort. 
Indeed a limit of the considered theoretical investiga- 
tion is that its systematic application to more complex 
cases is challenging. However, the set of techniques we 
employed can be successfully used to study a large 
class of simple systems, and are in general suitable for 
the analytical investigation of structural robustness of 
biological networks, complementary to simulations and 
experiments. 

Notes 

1 The concentration mismatch is more "softly" con- 
trolled, since the derivative of the difference 
x 2 — x x = d — C2 + b 22 x 2 — bnXi is not influenced by the 
nonlinear term bi 2 (xi, x 2 )x 2 = b 2 i(xi, x 2 )x\. 

2 I.e. the nullclines have no common tangent lines. 

3 Cf. the erratum: http://www.math.rutgers.edu/~son- 
tag/FTPDIR/angeli-ferrell-sontag-pnas04-errata. txt and 
[42]. 

Additional material 



Additional file 1: One additional file includes the proofs for 
Propositions 3, 4, 5 and 6 in the main paper. 



Acknowledgements 

The authors acknowledge financial support by the National Science 
Foundation (NSF) grant CCF-0832824 (The Molecular Programming Project). 
We are grateful to R. M. Murray, for helpful advise and discussions, and to 
the Reviewers for their constructive comments. 

Author details 

^ipartimento di Matematica ed Informatica, Universita degli Studi di Udine, 
Via delle Scienze 206, 33100 Udine, Italy. 2 Division of Engineering and 
Applied Science, California Institute of Technology, 1200 E. California Blvd. 
Pasadena, CA 91125, USA. 

Authors' contributions 

FB and EF performed research and wrote the paper. 

Received: 10 November 2010 Accepted: 17 May 2011 
Published: 17 May 2011 

References 

1. Kitano H: Systems biology: A brief overview. Science 2002, 
295(5560):1662-1664. 

2. Barkai N, Leibler S: Robustness in simple biochemical networks. Nature 
1997, 387(6636):913-917. 

3. Alon U, Surette MG, Barkai N, Leibler S: Robustness in bacterial 
chemotaxis. Nature 1999, 397(671 5):1 68-1 71. 

4. Yi TM, Huang Y, Simon Ml, Doyle J: Robust perfect adaptation in bacterial 
chemotaxis through integral feedback control. Proceedings Of The 
National Academy Of Sciences Of The United States Of America 2000, 
97(9):4649-4653. 

5. Ma W, Trusina A, El-Samad H, Lim WA, Tang C: Defining Network 
Topologies that Can Achieve Biochemical Adaptation. Cell 2009, 

138(4):760-773. 

6. Prill RJ, Iglesias PA, Levchenko A: Dynamic Properties of Network Motifs 
Contribute to Biological Network Organization. PLoS Biology 2005, 3(11): 
e343.. 



Blanchini and Franco BMC Systems Biology 201 1, 5:74 Page 14 of 14 

http://www.biomedcentral.eom/1752-0509/5/74 



10. 



12. 



13. 



14. 



16. 



17. 



19. 



20. 



22. 



23. 



24. 



25. 



26. 



27. 



28. 



29. 



30. 



32. 



33. 



Kwon YK, Cho KH: Quantitative analysis of robustness and fragility in 
biological networks based on feedback dynamics. Bioinformatics 2008, 
24(7):987-994. 

Gomez-Gardenes J, Y M, Fiona LM: On the robustness of complex 
heterogeneous gene expression networks. Biophysical Chemistry 2005, 
115:225-229. 

Gorban A, Radulescu O: Dynamical robustness of biological networks 
with hierarchical distribution of time scales. IET Systems Biology 2007, 
1(4):238-246. 

Kartal O, Ebenhoh O: Ground State Robustness as an Evolutionary Design 

Principle in Signaling Networks. PLoS ONE 2009, 4(12):e800L 

Aldana M, Cluzel P: A natural class of robust networks. Proceedings of the 

National Academy of Sciences of the United States of America 2003, 

1 00(1 5):871 0-871 4. 

Tian T: Robustness of mathematical models for biological systems. 

ANZIAM J 2004, 45:C565-C577. 

Shinar G, Milo R, Rodriguez Martinez M, Alon U: Input-output robustness 
in simple bacterial signaling systems. Proceedings of the National Academy 
of Sciences 2007, 104:19931-199935. 

Shinar G, Feinberg M: Structural Sources of Robustness in Biochemical 

Reaction Networks. Science 2010, 327(5971 ):1 389-1 391. 

Feinberg M: Chemical reaction network structure and the stability of 

complex isothermal reactors - I. The deficiency zero and deficiency one 

theorems. Chemical Engineering Science 1987, 42:2229-2268. 

Sontag E: Monotone and near-monotone biochemical networks. Systems 

and Synthetic Biology 2007, 1:59-87. 

Angeli D, Ferrell JE, Sontag ED: Detection of multistability, bifurcations, 
and hysteresis in a large class of biological positive-feedback systems. 

Proceedings of the National Academy of Sciences of the United States of 
America 2004, 1 01 (7):1 822-1 827. 

Jacquez J, Simon C: Qualitative Theory of Compartmental Systems. SIAM 
Rev 1993, 35:43-79. 

Abate A, Tiwari A, Sastry S: Box Invariance for biologically-inspired 
dynamical systems. 46th IEEE Conference on Decision and Control, New 
Orleans, LA 2007, 5162-5167. 

El-Samad H, Prajna S, Papachristodoulou A, Doyle J, Khammash M: 
Advanced Methods and Algorithms for Biological Networks Analysis. 

Proceedings of the IEEE 2006, 94(4):832-853. 

Radde N, Bar N, Banaji M: Graphical methods for analysing feedback in 

biological networks - A survey. Int J Syst Sci 2010, 41:35-46. 

Rouche N, Habets P, Laloy M: In Stability theory by Liapunov's direct method. 

Volume 22. New York: Springer-Verlag; 1977, [Applied Mathematical 

Sciences, xii+396 pp. ISBN 0-387-90258-9]. 

Blanchini F: Set invariance in control - a survey. Automatica 1999, 

35(1 1):1 747-1 767. 

Blanchini F, Miani S: In Set-theoretic methods in control. Volume 22. Boston: 
Birkhauser; 2008, [Systems & Control: Foundations & Applications]. 
De Jong H: Modeling and simulation of genetic regulatory systems: a 
literature review. Journal of Computational Biology 2002, 9:67-103. 
Alon U: An Introduction to Systems Biology: Design Principles of Biological 
Circuits Chapman & Hall/CRC; 2006. 

Nikolov S, Yankulova E, Wolkenhauer O, Petrov V: Principal difference 
between stability and structural stability (robustness) as used in systems 
biology. Nonlinear Dynamics Psychol Life Sci 2007, 1 1 (4):41 3-433. 
Mangan S, Zaslaver A, Alon U: The coherent feedforward loop serves as a 
sign-sensitive delay element in transcription networks. Journal of 
Molecular Biology 2003, 334:197-204. 

Levine E, Zhang Z, Kuhlman T, Hwa T: Quantitative Characteristics of Gene 

Regulation by Small RNA. PLoS Biology 2007, 5(9):e229.. 

Mitarai N, J B, S K, S S, Z C, E M, Sneppen K: Dynamic features of gene 

expression control by small regulatory RNAs. Proceedings of the National 

Academy of Sciences of the United States of America 2009, 

106(26):1 0655-1 0659. 

Mehta P, Goyal S, Wingreen NS: A quantitative comparison of sRNA-based 
and protein-based gene regulation. Mol Syst Biol 2008, 4. 
Williamson T, Schwartz JM, Kell DB, Stateva L: Deterministic mathematical 
models of the cAMP pathway in Saccharomyces cerevisiae. BMC Systems 
Biology 2009, 3. 

Jacob F, Perrin D, Sanchez C, Monod : L'opeeron: groupe de genes a 
expression coordonnee par un operateur. J C R Acad Sci 1960, 
250:1727-1729. 



34. 



35. 



36. 



37. 



38. 



39. 



40. 



42. 



Vilar JMG, Guet C, Leibler S: Modeling network dynamics: the lac operon, 
a case study. Journal of Cell Biology 2003, 161(3):471-476. 
Yildirim N, Mackey M: Feedback Regulation in the Lactose Operon: A 
Mathematical Modeling Study and Comparison with Experimental Data. 

Biophysical Journal 2003, 84(5):2841-2851. 

Yildirim N, Santillan M, Horike D, Mackey M: Dynamics and bistability in a 
reduced model of the lac operon. Chaos 2004, 14(2):279-292. 
Huang CYF, Ferrell JJ: Ultrasensitivity in the mitogen-activated protein 
kinase cascade. Proceedings Of The National Academy Of Sciences Of The 
United States Of America 1996, 93:10078-10083. 

Ferrell J, James E, Machleder EM: The Biochemical Basis of an AII-or-None 
Cell Fate Switch in Xenopus Oocytes. Science 1998, 280(5365):895-898. 
Markevich Nl, Hoek JB, Kholodenko BN: Signaling switches and bistability 
arising from multisite phosphorylation in protein kinase cascades. J Cell 
Biol 2004, 164(3):353-359. 

Qiao L, Nachbar RB, Kevrekidis IG, Shvartsman SY: Bistability and 
Oscillations in the Huang-Ferrell Model of MAPK Signaling. PLoS Comput 
Biol 2007, 3(9):e184. 

Ferrell JE, Pomerening JR, Kim SY, Trunnell NB, Xiong W, Huang CYF, 
Machleder EM: Simple, realistic models of complex biological processes: 
Positive feedback and bistability in a cell fate switch and a cell cycle 
oscillator. FEBS letters 2009, 583(24):3999-4005. 

Russo C, Giuraniuc C, Blossey R, Bodart JF: On the equilibria of the MAPK 
cascade: Cooperativity, modularity and bistability. Physica A: Statistical 
Mechanics and its Applications 2009, 388(24):5070-5080. 



doi:1 0.1 1 86/1 752-0509-5-74 

Cite this article as: Blanchini and Franco: Structurally robust biological 
networks. BMC Systems Biology 201 1 5:74. 



Submit your next manuscript to BioMed Central 
and take full advantage of: 

• Convenient online submission 

• Thorough peer review 

• No space constraints or color figure charges 

• Immediate publication on acceptance 

• Inclusion in PubMed, CAS, Scopus and Google Scholar 

• Research which is freely available for redistribution 



Submit your manuscript at 
www.biomedcentral.com/submit 



o 



BioMed Central 



